Mon. Not. R. Astron. Soc. OOP. [TlfTOl C2002') Printed 2 February 2008 (MN I*TeX style file v2.2) 



Long-term Nonlinear Behaviour of the Magnetorotational 
Instability in a Localised Model of an Accretion Disc 

L. J. Silvers* 

Laboratoire de Radioastronomie, Departement de Physique Ecole Normale Superieure, 24 Rue Lhmond, 75231, Paris, 05, Prance. 



ABSTRACT 

For more than a decade, the so-caUed shearing box model has been used to study 
the fundamental local dynamics of accretion discs. This approach has proved to be 
very useful because it allows high resolution and long term studies to be carried out, 
studies that would not be possible for a global disc. 

Localised disc studies have largely focused on examining the rate of enhanced 
transport of angular momentum, essentially a sum of the Reynolds and Maxwell 
stresses. The dominant radial-azimuthal component of this stress tensor is, in the 
classic Shakura-Sunayaev model, expressed as a constant a times the pressure. Pre- 
vious studies have estimated a based on a modest number of orbital times. Here we 
use much longer baselines, and perform a cumulative average for a. Great care must 
be exercised when trying to extract numerical a values from simulations: dissipation 
scales, computational box aspect ratio, and even numerical algorithms all affect the 
result. This study suggests that estimating a becomes more, not less, difficult as com- 
putational power increases. 
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1 INTRODUCTION 

I Linear stability analysis has shown that a differentially ro- 
' tating disc, with angular velocity decreasing outwards, is 
, unstab le in the presence o f a weak magnetic field (se e , for e x- 
■ ample. iBalbus fc HawlevI l|l99ll ): iBalbus fc HawlevI l(l992h ). 
' This magnetorotational instability (MRI) gives rise to tur- 
bulence within the disc and is a plausible way that efficient 
outward transport of mo mentum in a disc can be realised 
l|Balbus fc HawlevI l|l998l )V However, while this is an appeal- 
ing concept, it is crucially important, due to the non-linear 
form of the governing equations, to investigate the non-linear 
evolution and saturation of this instability to determine if 
this is indeed a route by with angular momentum transport 
can be achieved on the required time-scales. 

The equations that govern the evolution within a disc 
are computationally demanding to evolve. Therefore, in or- 
der to study the evolution and saturation of key quantities 
associated with this instability (e.g. stresses, energies etc), 
it is useful to e mploy a localised, 'shearing-box' model (see, 
for e x ample, iBrandenburg. Nor dlund. S tein fc TorkelssonI 
j 19951 '): iHawlev. Gammie fc Balbus (1995, )). This approach 
has proved useful as it has enabled higher local resolutions to 
be achieved for the same computational cost as that for a full 
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disc simulation. However, while a localised approach does fa- 
cilitate longer runs and higher resolutions, at the point of 
the early shearing box calculations the resources were such 
that only modest run times (u sually less than 30 orbits) 
were p ossible (see, for example Hawley, Gammie fc Bal bua 

Hiil)). 

When discussing the MRI, the principal quantity of 
interest is a, which is related to the rate of angular 
momentum trans port in a di s c. This quantity is defined 
as (jShakura fc SunvaevI l| 19731 ): iHawlev. Gammie fc Balbu3 
119951)): 

a = Wxy/P (1) 

In this definition, P is the gas pressure and Wxy is a compo- 
nent of the X — y stress tensor i.e. 

BxBy , , 

Wxy = pVxVy , (2) 

47r 

where p denotes density; Vx is the x component of the veloc- 
ity field; Vy is the perturbed component of the velocity field 
in the y direction; Bx is the x component of the magnetic 
field and By is the y component of the magnetic field. Q 
In turbulent flows transport quantities such as a are 

^ Note here that (x,y) refers to a local set of coordinates in (r, </>) 
and not [—(j>,r) as sometimes utilised in certain hydrodynamical 
models. 
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highly fluctuating. The resource constraints at the time of 
the earliest calculations forced a values to be quotes from by 
a simple average over a modest number of orbits, which is 
useful to give an idea o f the order of magnitude of this qua n- 
tity (see, for example. iHawlev. Gammie fc Balbui (1199^) ). 
However, unless an extremely large number of orbits is con- 
sidered, the saturation level will vary as you increase the 
number of orbits over which the average is evaluated. This 
can, make it difficult to see clearly the effect if varying each 
of the parameters (e.g. dissipation scale, domain size etc) as- 
sociated with the problem. It is better, if the resources are 
available, to calculate a cumulative average for this fluctu- 
ating quantity. This does require vast numbers of computer 
hours, which is only recently possible. As such, the princi- 
pal aim of this paper is to calculate long-time cumulative 
averages of a. 



In order to keep the computational costs acceptable, 
we choose to make several simplifying assumptions in this 
paper. First, we choose to initially work in a cubic do- 
main, which is only relaxed later in the paper. This en- 
abled us to explore the effect of decreasing the dissipa- 
tion scale on the cumulative average a and to show that 
small scales are playing a role in determining the satu- 
ration level of this quantity. The second simplification is 
that we chose to solve the ideal MHD equations. This re- 
duces the cost of the calculation by a significant fraction 
though it does mean that all calculations are at an effec- 
tive Prandtl number of order unity. This is acceptable as 
in this particular investigation we were not trying to ex- 
amine the effect of varying the Prandtl number. I ndeed , 
this has recently been considered llLesur fc Logarettil l|2007h : 
iFromang. Papaloizou. Lesur. HeinemannI 1 20081 )'). The work 
in this paper is aimed, somewhat, to complement these 
recent works even though here we do not explicitly in- 
clude dissipation c oefficients. We will show that, as in the 
iLesur fc Logarettil (2007) paper the effect of increasing reso- 
lution for a fixed box size is to increase a. Further, we show 
that the effect of increasing the box length is to decrease the 
saturation level of a as the parasitic instability is allowed to 
take hold. 



The paper will proceed as follows: In section 2 we will 
fully detail our model and the numerical method used in 
the solution of the governing equations. We will discuss the 
results in section 3. This section is divided into two subsec- 
tions. In section 3.1 we begin by showing the problems that 
are encountered in standard shearing box calculations that 
makes long-term cumulative averages difficult to calculate. 
We show that this problem is a result of compounded errors 
that come from the shearing boundary conditions. We then 
present a correction to the numerical algorithm, which per- 
mits us to calculate long-term cumulative averages and show 
how Q varies with resolution (or dissipation length scale) . In 
section 3.2 we bring the results section to a close by show- 
ing that the effect of increasing the length of the box in 
the azimuthal direction is to decrease a. In Section 5 we 
summarise the finding as well as discuss the astrophysical 
consequences. 



2 THE MODEL 

In this paper we restrict consideration to a local patch of a 
Keplerian disc and to use the shearing box approximation. 
Within this paper we choose to follow the strategy of many 
of the earlier works by which we mean that we will consider 
an ideal MHD problem. We shall assume for simplicity in 
this paper that the z-component of gravity is suitably small 
and so may be neglected (assumed to be within a pressure 
scale height of the midplane of the disc) . Further we assume 
that the gas is isothermal. As such, equation set may be 
written as: 



dp 
dt 



+ V.pv = 



(3) 



'dt 



+ v.Vv = 



V P + 
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Stt 
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V.B = 



P = 



(4) 



(5) 



(6) 



(7) 



where we obey the standard notational convection and thus, 
p denotes density, P denotes the gas pressure, v and B are 
the velocity and magnetic fields, respectively, Q. is the angu- 
lar velocity, q = —d\nQ./d\nR. 

Our work will differ from earlier works in the follow- 
ing ways. First we choose to start from a perturbation that 
has a specific analytic form that is not a random perturba- 
tion at each of the grid points. This initial perturbation in 
such simulations are immaterial so long as they are small in 
amplitude. However, in order to be completely clear on our 
starting condition, we choose to start from a perturbation 
that can be expressed in a simple analytic manner so that 
we have a consistent start condition as we vary the resolu- 
tion. We choose to perturb each of the velocity components 
in the following manner: 

Defining 

5 

= 5Z sinpi-rix - (1.2 - 0.2n)], (8) 



fy = sin[2™y - (1 - 0.2n)], (9) 



f, = J2 sin[27rn2 - (0.8 - 0.2n)] (10) 

11=1 

the initial velocity components are then expressible as: 

= o.omf^fyf, (11) 



Vy = -qnX + O.Oinf.,fyf, 



(12) 



V, = omnf^fyf, (13) 

We have chosen this particular form to give the perturbation 
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Table 1 . Summary of all the cases presented in the paper 



Case 


Resolution 


Aspect 


Artificial 


Mean 










Ratio 


Viscosity 


Cleaning 


1 


32 


32 : 


32 


1:1:1 


No 


No 


2 


64 


64 : 


64 


1:1:1 


No 


No 


3 


128 : 


128 


: 128 


1:1:1 


No 


No 


4 


32 


32 : 


32 


1:1:1 


Yes 


No 


5 


64 


64 : 


64 


1:1:1 


Yes 


No 


6 


32 


32 : 


32 


1:1:1 


No 


Yes 


7 


64 


. 64 : 


64 


1:1:1 


No 


Yes 


8 


128 : 


128 


: 128 


1:1:1 


No 


Yes 


9 


64 : 


128 


: 64 


1:2:1 


No 


Yes 


10 


64 : 


192 


: 64 


1:3:1 


No 


Yes 


11 


64 : 


384 


: 64 


1:6:1 


No 


Yes 



some complexity while still keeping it simple to code and 
express. 

In this work we choose to start from a uniform density 
and so apply no perturbation to this quantity. We take po = 
1, Po = 10~^, n = Cs = 10~^ and q is taken to be that for a 
Keplerian profile namely, 3/2. The embedded magnetic field 
is uniform, in the vertical-direction and such that /3 = 800. 

Initially, unlike in most of the earlier works we choose 
to consider a shearing box where each side is of unit length 
i.e. one that has the same length and number of grid points 
in each directions. Such an approach facilitates probing to 
higher resolutions without resorting to earlier methods of 
assuming apriori that less resolution is required in the y- 
direction. While less resolution might be warranted in the 
y-direction in simulations that have a clear difference in the 
scale of the structures in each of the spacial directions we 
choose here to use caution and not make this assumption 
which, amongst other things, gives rise to an anisotropy in 
the numerical dissipation. 

The equations, given above, are solved with a parallel 
finite-difference time-explicit code. The code is essentially 
just a parallelized version of the serial Z EUS code that is 
freely available. As in earlier papers (e.g. IStone fc Norman 
'l992 [) : lHawlev fc Stond (| 19951 ) : lHawlev. Gammie &: Balbus 



1996f )l we use the method of constrained transport scheme 
for magnetic advancement so that the V.B = is enforced 
to the accuracy of the machine. In this work we use the stan- 
dard shearing-box boundary conditions. Thus the domain is 
periodic i n both the y and z directions and s hearing-periodic 
in X (see iHawlev. Gammie fc Balbuj (Il996l ) for further de- 
tails) .This code contains an adaptive step size that is multi- 
plied, as is standard, by a safety factor. We have conducted 
several tests to ensure that the step size does not give rise to 
numerical instability. This check involved varying the safety 
factor in this code to ensure that all prominent global fea- 
tures, trends and saturation levels are fully captured by the 
integration, which could be referred to as a check that we 
have the 'true' or 'converged' solution. 



3 RESULTS 

Before beginning to discuss the finding it helpful to provide 
an overview of all cases at this point. Table [1] shows each of 
the cases that we have considered. We have used the method 



detailed in iHawlev. Gammie fc BalbusI l| 19951 ) to ensure that 
the fastest growing wavelength is resolved. 

We will choose to split the results section into two 
halves. Firstly we will discuss the effect of increasing the 
resolution and then we will go on to discuss the effect of in- 
creasing the box length, in the y direction, while maintaining 
the same resolution per unit length. 

3.1 Ideal MHD 

We begin the discussion of the results with our fiducial 
model, case 1. In this model we consider a resolution of 
32 grid points in each of the three directions. This fiducial 
model is at a similar resolution, in the x and y directions, as 
can be found in many of the original calculations. This said, 
the earlier calculations did typically a longer box length than 
we do here. This deliberate change was made for two rea- 
sons. First, the decrease in the box size enables us to push 
to much higher resolutions. Second, and most importantly 
for what we wish to examine here, an increased integration 
time i.e. such a reduction in the box length facilitates being 
able to ascertain many more orbits for the same computa- 
tional cost. We shall discuss the ramifications of increasing 
the box length later in the paper and some comments on the 
effect of varying the compu t ation al domain can be found in 
[Winters. Balbus fc Hawlevl (|2003l ). 

As stated in the introduction, we wish to determine a 
cumulative average for the quantity a, which is related to 
the difference of the Reynolds and Maxwell stresses. Such a 
quantity can only be calculated if the simulations are run for 
sufficiently long so that there is a large time in the saturated 
state for the cumulative average to settle. Figure [1] shows the 
initial section evolution of the magnetic and kinetic energies 
for the 32^ case and illustrates the it can take several orbits 
for the saturated state to begin. We note that, despite the 
reduction in the box size, these figures resemble those from 
the earlier calculations. However, we note that time to reach 
a statistically steady state occurs at a later point in terms 
of the number of orbits but at a comparable time from the 
start of the simulations. Other slight differences from earlier 
calculations are combined effect of box size, isothermality 
and choice of plasma (3. 

These figures appear to show that the evolution is com- 
plete and that the final state has been obtained and thus one 
might thing that all that remains to calculate a cumulative 
average is to carry the calculation on for sufficient time and 
then evaluate the cumulative average of a. However, figure 
[2] shows that, in this calculation, what appears to initially 
be a statistically steady state is only approximately so. If 
you examine the solution for 1000 orbits one sees that the 
first 200 orbits or so is in fact a slow ramp up to a higher 
saturation level. 

The most important implication from a shift in the av- 
erage magnetic and kinetic energies is that this is almost 
certainly a corresponding shift in the momentum transport 
and figure [2] shows that this is indeed the case. However, 
from a figure of this kind it is exceptionally hard to gain a 
feeling for the structure — indeed much of the information 
on the temporal variability is being lost. As such, we provide 
a more detailed series of snapshots over the full 1000 orbits 
in |3l This images show that the system is moving from a 
state which is rather turbulent into a state where there is a 
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Figure 3. A series of plots sliowing in more detail the evolution of a as the number of Orbits increases for case 1 
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Figure 1. Top: Kinetic energy evolution for tlie fiducial case 
for the first 50 Orbits. Bottom: Magnetic energy evolution (nor- 
malised) for the fiducial case for the first 50 orbits. 



repeated classical channel solution presenting before being 
destroyed and the cycle repeating. 

In the light of the images show above, one of the most 
natural questions concerns whether the same behaviour is 
still found if the length-scale for the dissipation scale is de- 
creased. In this problem this is akin to increasing the number 
of grid points that are used. As such we have examined the 
case when the number of grid points in each spatial direction 
and the findings are shown in figure [l] Note that, while the 
chmb to the saturated state is slightly slower in this case, 
a transition from one quasi-statistically steady state to an- 
other still occurs. Similar was found at a resolution of 128^ 
and so we conclude that this behaviour is independent of 
the length-scale for dissipation. 
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Figure 2. Top: Kinetic energy evolution (normalised) for the case 
1 for 1000 orbits. Middle: Magnetic energy evolution (normalised) 
for case 1 for 1000 orbits. Bottom: The evolution of a for case 1. 



3.1.1 Artificial Viscosity 

In the above calculations we have not included a resistivity 
or a viscosity of any kind and have just left the dissipation 
to be the result of the grid. However, there in several of 
the local shearing box calculations the was an artificial vis- 
cosity term included in the e quati on of motion of the form 
l|Hawlev. Gammie fc BalbusI l|l996l )): 



A^Cp{S,v,){Ax) (14) 

The motivation for such a term was to provide lim- 
ited extra dissipation in the case when there is a strong 
compressional solution. As such, it retains the original ideal 
magnetohydrodynamic (MHD) nature of the equations and 
only calls for extra assistance is special situations. While we 
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Figure 4. Top: Kinetic energy evolution (normalised) for case 2 
for 1000 orbits. Middle: Magnetic energy evolution (normalised) 
for case 2 for 1000 orbits. Bottom:The evolution of a for the case 
where there are 64'^ grid points (case 2). 
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Figure 5. The evolution of a for the case where there are 32^ 
grid points and an artificial viscosity (case 4). 



about might in some way be the result of the lack of the 
inclusion of an artificial viscosity. As such, for completeness 
we here examine the 32'' calculation with the artificial vis- 
cosity term included and where C — 2.0 as is standard. 
This implies that smoothing is carried out over two cells. 
F urther dis cussion of the artificial viscos ity can be found 



Balbus (199^) and IStone fc NormanI 



__.Hawlev. Gammic 
(1 19921) . 

Figure [S] shows that the same global trend is witnessed 
in this case as in the case where there is no artificial viscos- 
ity included. This shows that the inclusion of an artificial 
viscosity does not change the result found above. 



3.1.2 Problems with the Shearing Boundary Conditions 

The behaviour witnessed above is fully explained by an ex- 
amination of what occurs to the mean of the various com- 
ponents. Many of these components should formally be zero 
but the interpolation necessary in the implementation of the 
shearing boundary conditions give rise to errors as shown 
in figure |6] Thes e errors have been mentioned be fore in 
earlier calculations iHawlev. Gammie fc BalbusI (|l995l ). How- 
ever, many of the early calculations were for short time spans 
and there the error was stated to be about 10~^Cs. while we 
confirm that this is true for small numbers of orbits this is 
not true in the long term. In fact, as figure|6]shows the mean 
Bz and By fields both grow and that the By field grows to 
such and extent that it becomes of greater magnitude. We 
note here that the numerical method employed prevents Bx 
from obtaining significant error. The constrained transport 
forces Bx to stay constant to round off because it is evolved 
with the emf components Ey and which are strictly peri- 
odic and hence Bx is conserved to round off error. However, 



would argue that a viscosity of some kind is useful we would 
argue that it is much better to consider a physical viscosity 
that acts continually. 

However, while we do not believe that the artificial vis- 
cosity method is the most instructive it is a method that 
has served the community for many years and, as such, 
one must entertain the possibility that the solutions found 



interpolation in Ex. 

If we are to obtain a valid long-term average for a then 
we need to remove this error, which compounds to give us 
an unintentional switch in the mean field. As such we pro- 
ceed by forcing all mean values to their mathematically for- 
mal values by removing the error at set intervals. Obviously 
the most rigorous correction would be invoke this cleaning 
procedure every time-step. However, this is computationally 
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Figure 6. The evolution of each of the components of the mag- 
netic field for case 1. 



Figure 9. The evolution of a when full cleaning is employed for 
case 7 (1:1:1 aspect ratio) 





1.8 




1.6 


o 




CL 


1.4 














1.2 


CQ 




m'^ 


1 


V 




1 

A 


0.8 


>^ 




>^ 


0.6 








0.4 




0.2 








data 
mean 
-running average 




i-^-i s- It;: t" fe- W <i 



200 



400 600 
Orbits 



800 



1000 



Figure 7. The evolution of a when full cleaning is employed for 
case 6 (1:1:1 aspect ratio) 



Table 2. Selected cases together and the associated a from a 
cumulative average 



Case 


Resolution 


Box dimension 


a 


6 


32:32:32 


1:1:1 


0.85 


7 


64:64:64 


1:1:1 


0.94 


8 


128:128:128 


1:1:1 


1.02 


9 


64:128:64 


1:2:1 


0.43 


10 


64:192:64 


1:3:1 


0.29 


11 


64:384:64 


1:6:1 


0.16 



in table [2] (cases 6-8) , which shows a clear upward trend as 
resolution is increased. As such, this shows that small scales 
are important in these calculations and the level at which 
the dissipation scale is set is crucial is determining that rate 
of enhanced transport. 



demanding. Therefore, we choose here to correct the mean 
values every 100 time-steps as we determined that this gives 
phenomenally good agreement with cleaning every time-step 
for considerably less cost. 

Figure [7] shows the new temporal evolution of a once 
the cleaning has been carried out for the 32'^, which is case 
6. It is clear that the effect of cleaning the mean values on 
a is to continue with the saturated level that we witnessed 
in the first part of figure [2] Therefore, it is now possible to 
evaluate the cumulative average until it saturates. Figure [7] 
not only shows the cumulative average but also the mean 
value for this case. Figure |8] shows a more detailed series 
of snapshots for a in this case and clearly shows that the 
randomised, turbulent, behaviour persists for all time and 
we do not now move to a repeated channel solution case. 

The primary aim of this paper was to examine the effect 
of decreasing the dissipation scale (by increasing the resolu- 
tion) and we are now in a position to do so. Figure [9] show 
the effect of doubling the resolution (decreasing the dissipa- 
tion scale) on a. We note that the average value is increased 
as resolution is increased. This fact is amplified by the data 



3.2 The Effect of Increasing the Box length 

In the above section we chose to focus on a domain where 
the aspect ratio is 1 : 1 : 1. This was selected to facilitate 
the longer time calculations at higher resolutions. However, 
as the earlier calculations considered a 1 : 27r : 1 and so in 
this section we examine what happens at modest resolutions 
to the saturated level of a. We have considered boxes with 
the aspect ratios 1 : 2 : 1, 1 : 3 : 1 and 1:6:1. The effect 
of elongating the is summarised in [2] (cases 7, 9-11) and the 
evolution is shown in figure [TOl 

There is a clear decrease in the saturation level as the 
box is increased in the azimuthal direction. What we are 
seeing is the effect that i n short box leng t hs th e parasitic 
instability is limited (see iGoodman fc Xu (|l994l ) for more 
details on this instability). 



4 CONCLUSIONS 

In this paper we have returned to the original concept to 
understand turbulent transport of angular momentum in an 
accretion disc. As such, we considered the standard local 
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Figure 8. A series of figures showing the time evolution of a when 'cleaning' is employed for case 
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Figure 10. Top: The evolution of a when full cleaning is em- 
ployed for case 9 (1:2:1 aspect ratio). Middle:The evolution of a 
when full cleaning is employed for case 10 (1:3:1 aspect ratio). 
Bottom: The evolution of a when full cleaning is employed for 
case 11 (1:6:1 aspect ratio) 



shearing box model in which the ideal MHD equations were 
solved. We chose to start by considering a cube with sides 
of equal length because we wished to have large integra- 
tion times (1000 orbits) and a numerical dissipation that 
the same in all directions. 

Our calculations shown that, without correction, there 
is a shift in the accretion rate (a shift in the turbulent trans- 
port rate). Such a shift is interesting due to the fact that 
there are shifts in spectra of accretion discs, which imply a 



corresponding shift in accretion rate. However, in this case, 
the shift is unintentional and a by-product of the interpola- 
tion in the shearing box boundary conditions that introduce 
errors because of the interpolation that is a necessary part of 
the procedure. This error is unavoidable in such calculations 
and does remain small for low numbers of orbits. However, 
for large numbers of orbits the error in shown to compound 
to such an extent that the mean By field, which should be 
zero, grows and becomes larger in magnitude that the im- 
posed field. This error filers into the various components 
of both the velocity and magnetic fields. 

In this paper we sought to consider the simple case of an 
imposed vertical magnetic field i.e. without a shift to a state 
where where was a mean field in more than one direction. 
As such we applied a simple fix to this issue that removes 
the error that is induced into the mean value of each of 
the components of the velocity and magnetic fields. This 
facilitated us to evaluate a long-term cumulative average for 
Q, which is a better quantity when you wish to compare a 
values for different parameters. Following this procedure, we 
showed that increasing the resolution, for a fixed box size, 
gives rise to an increase in the value of a, which suggests that 
the values here are probably a lower bound for a values in 
these kind of simulations. This leads us to deduce that small 
scales are playing and important role in determining the rate 
of accretion. In an accretion disc the dissipative length-scale 
is much less than we considered here and so we conclude 
that in a disc there will be a significant enhancement to the 
rate of angular momentum transport. 

As mentioned above, the work in the first part of this 
paper is not in the standard domain of the early work in 
a shearing box i.e. it was conducted in a cube with sides 
of equal length. This was selected to facilitate our desire 
for higher resolutions and and extremely long integration 
times. However, it of obvious interest to examine the ef- 
fect of increasing the box length at modest resolutions. We 
found that increasing the box length led to a reduction in the 
mean value of a. This is attributed to the increased freedom 
for the parasitic instability as the box length is increased 
(see [Good man fc Xu {WM) for a detailed discussion of this 
instability). We note here that while that increasing the res- 
olution by a factor of four leads to a modest (15box length 
is much more profound and increasing the box length by a 
factor of three leads to a reduction of a such that it is almost 
2/3rds of its value for a box with length 1 in the y-direction. 

While we acknowledge that there is still considerable 
working to be conducted in this area to fully understand 
the magnetorotational instability we believe that this pa- 
per is an important contribution to the field as it clearly 
demonstrates the need for great care in numerical calcu- 
lations. Further it shows that, for the case of a imposed 
vertical magnetic field, a is non- negligible and it points to 
there existing sufficient transportation within a disc to give 
accretion in an appropriate time-scale. 

The reader should note that what has been found 
here is for the case of a uniform vertical magnetic 
field and that this is very different to the case of a 
zero-mean magnetic field as co n sidered by, for exam- 
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le, Pessah. Chan. Psaltis 1 ||2007|) Fromang fc Papaloizoul 



200i ) and lFromang. Papaloizou. Lesur. HeinemannI (|2008l ). 
where a in many situations was found to be extremely small. 
Here we have shown that a increases as the scale of dissipa- 
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tion decreases and that it decreases as the length of the com- 
putational domain in the azimuthal direction is increased. It 
is not clear to what extend zero-mean flux calculations for 
long times are effected by the error that we found and we 
would suggest that all calculations in a shearing box check 
that mean quantities are maintained to expected values. 

This paper does pose some questions of which the most 
obvious is whether or not a change in the balance of moan 
components of the magnetic field, along with the strength 
of the magnetic field, always gives rise to a change in the 
accretion rates. This issue is being investigated at present. 
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